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ABSTRACT 



Context. The early phase of the coalescence of supermassive black hole (SMBH) binaries from their host galaxies provides a guar- 
anteed source of low-frequency (nHz-/jHz) gravitational wave (GW) radiation by pulsar timing observations. These types of GW 
sources would survive the coalescing and be potentially identifiable. 

Aims. We aim to provide an outline of a new method for detecting GW radiation from individual SMBH systems based on the Sloan 
Digital Sky Survey (SDSS) observational results, which can be verified by future observations. 

Methods. Combining the sensitivity of the international Pulsar Timing Array (PTA) and the Square Kilometer Array (SKA) detectors, 
we used a binary population synthesis (BPS) approach to determine GW radiation from close galaxy pairs under the assumption that 
SMBHs formed at the core of merged galaxies. We also performed second post-Newtonian approximation methods to estimate the 
variation of the strain amplitude with time. 

Results. We find that the value of the strain amplitude h varies from about 10~ 14 to 10~ 17 using the observations of 20 years, and we 
estimate that about 100 SMBH sources can be detected with the SKA detector. 

Key words. Binaries:general -Black hole physics -Stars: evolution - gravitational waves 
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A gravitational wave (GW), which is described as a space per- 
turbation of the metric traveling at the speed of light, is a natural 
consequence of Einstein's theory of gravity ( 'general relativity' , 
Einstein 1916, 1918). It has been accepted (Thorne & Braginskii 
1976) that supermassive black hole (SMBH) binaries, which 
are ubiquitous in the nuclei of low-redshift galaxies (Magorrian 
et al. 1998), are boua fide GW sources, with the international 
Pulsar Timing Array (PTA) and the proposed Square Kilometer 
Array (SKA) (Detweiler 1979; Hellings & Downs 1983; Kaspi 
et al. 1994; Wyithe & Loeb 2003; Jenet et al. 2005; Yardley 
et al. 2010). The precision of the rotation periods for millisec- 
ond period pulsars, which is established via pulse arrival time 
measurements, allows the detection of stochastic SMBH GW 
background spectrum at nHz frequencies (Rajagopal & Romani 
1995; Jenet et al. 2006). These observations give us a good 
prospect for detecting GWs from the ensemble of these SMBH 
coalescence events throughout the universe. Some of these stud- 
ies have been taken up by other groups over the years. In partic- 
ular, the sensitivity curve of the SKA to GWs emitted by SMBH 
systems has been obtained from 100 pulsars using the Australia 
Telescope National Facility pulsar catalog (Yardley et al. 2010). 
The frequent growth of galaxies by merger processes and 
SMBH in galaxies result in some obvious mechanisms about 
the formation of SMBHs. Several studies have argued that an 
occupation fraction of SMBH as low as 0.01 at reshift 5 can 
be obtained from subsequent mergers (Menou et al. 2001 and 
references therein). Spherical galaxies often lead to coalescence 
timescales that are longer than the Hubble time, while highly 



flattened or triaxial galaxies lead to faster coalescence (Wyithe & 
Loeb 2003). Therefore this expected GW radiation investigation 
from coalescing SMBHs depends on the merger rate of massive 
galaxies, the demographics of SMBHs at low- and high- redshift, 
and the dynamics of SMBH binaries. Meanwhile the SMBHs 
need to reach a regime where GW radiation can drive the main 
evolutionary process (Begelman et al. 1980). If the galaxy pairs 
can evolve to this GW radiation regime, they will gradually lose 
angular momentum to GW radiation and will eventually coales- 
cence. The emission of GWs from SMBH sources requires that 
the binary coalesce in less than the Hubble time. 

Some results of close galaxy pairs have been released from 
the Sloan Digital Sky Survey (SDSS) data (Ellison et al. 2008, 
2010), such as the star formation rate, the merger rate of massive 
galaxies, the demographics of SMBHs at low redshift (z < 0.1) 
and the dynamics of SMBHs. A binary population synthesis 
(BPS) approach has been applied to study the characteristics of 
clusters and galaxies (see, e.g., Han et al. 2007; Zhang et al. 
2010). This method has also been systematically taken up by 
some groups to investigate the GW radiation from compact bi- 
naries in the Galaxy (Nelemans et al. 2001, Liu 2009; Liu et 
al. 2010; Ruiter et al. 2010). However, the GW radiation from 
SMBHs, combined with the BPS model and observation results, 
are not well known. Here we report how BPS , using SDSS re- 
sults, can be used to determine the GW radiation from SMBHs. 

In the next section, the BPS model is described. In Sect. 3, 
we present the results of our simulations. And the discussion is 
given in Sect. 4. 
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2. The binary population synthesis model 

2.1. Confirming the evolution trace of SMBH systems with 
OJ287 

Supermassive black holes are expected to be a type of GW 
sources, but few studies can accurately determine a galaxy 
pair's internal compositions and kinematics equations. The BL 
Lacertae object OJ287 (Sillanpaa et al. 1988, 1996; Valtonen 
2007) provides a possible insight into the precessional elliptical 
orbital evolution through GW radiation in the phase of coales- 
cence of SMBH binaries. The basic elements of this phase are as 
follows: 12-year intervals in optical outburst, component masses 
1.3 x 10 8 M o and 1.8 x 10 10 M G , an orbital period of nine years, 
an eccentricity of 0.67, a distance of 1.3 Gpc. 

In accordance with Einstein's general relativity, the metric 
can be written as 
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where r\ a p is Minkowsky metric, and h a p stands for GWs with 
I h a p |<s 1. The amplitude and the luminosity of GW radiation 
are 
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where G and c are the gravitational constant and speed of light, 
respectively. % is the distance from the GW sources to the Earth, 
Dij is the polar moment of mass, and the symbol <> denotes the 
averaging. 

The SMBHs more increasingly closer because of GW radi- 
ation over time, therefore the lifetime of galaxy pairs before the 
final in-spiral is (Jaffe & Backer 2003) 



fGW = 1.1 xl0 6 yrM, 
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where Mg stands for the total mass of the binary in units of 
10 8 M G , q (= mi/m.2 < 1) is the mass ratio, P or b is the orbital 
period in years. Meanwhile, the orbital angular momentum loss 
from two point masses is given by 
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where r is the separation of SMBHs in the galaxy pairs. 

As displayed above, an important quantity to be considered 
is the GW amplitude h, here we give a simple description of the 
calculation of h in our simulations. First, we give the definition 
for the polar moment of mass D t j in Eq. (3) 



L>ij - qij - -OijO qui, 



(6) 



where qij - 3 Jy'y^T°°(t,y), and i, j = 1,2. In general, we need 
to calculate this energy momentum tensor for an isolated system. 
For our particle binaries we just simulate the equation of motion 
of each BH using a numerical method. The primary of an SMBH 
system is defined here as the initially more massive component 
(mi) throughout. We also assumed that the primary BH m 2 is at 
the origin of coordinate, and that the secondary BH m\{< mi) 
traces out the precessional elliptical orbit around this primary. 



Fig. 1. Orbital trace template of SMBH systems using OJ287. 
The primary BH is at coordinates (0,0). The precessional ellipti- 
cal curve represents the evolutionary trace of the secondary BH. 



According to Eq. (16) of Antonacopoulos (1979), the equa- 
tion of motion of the secondary m\ can be described as 
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where x i and a: 2 are the vector position of secondary and primary 
in the BH binaries, respectively. And u \ is the vector velocity of 
secondary BH. To facilitate the calculation, we derive from Eq. 
(7) using the polar coordinate 
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here r = \x\ — x 2 \. Here we neglect the influence of the smaller 
quantity (such as 0(c 5 )) in Eqs. (7)-(9). We did not perform the 
calculation in full relativity. 

To fix the orbital trace of SMBHs, which is consistent with 
OJ 278's modulation, we investigated the numerical result us- 
ing the followed initial parameters of the position of periastron: 
ro = 3162AU, <po - 0, ro = and the tangent components of the 
velocity vo = 0.275c, the subscript denotes the initial values. 
Therefore, Fig. 1 shows an orbital evolution trace of the sec- 
ondary BH. Note that our calculation of the precessing-binary 
model is slightly different from that of Fig. 1 of Valtonen (2007), 
which is consistent with the fixed optical outbursts of OJ287 ac- 
cording to an accretion disk model. In other words, we excluded 
the influence of the accretion disk from our definition. That is 
mainly because the accretion disk dominates the electromagnetic 
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radiation. Our simulations focus on the importance of GW radi- 
ation. 

To analyze the GW amplitude over time, from Eq. (6), the 
components of the mass tensor (Capozziello et al. 2008) are 



D n = mir 2 (3cos 2 ^ - 1) 
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Here other quantities Djj - 0. The change of the GW amplitude 
is illustrated in Fig. 2. We can note four kinds of information 
from it. Firstly, we consider OJ287 as our standard template and 
obtain some typical feedback: the orbital period is 9.4 yr, the pre- 
cessional motion is 40° per period, and the eccentricity is 0.71. 
Secondly, this numerical curve covers about 40 years from 1996 
to 2036. And the maximum quantities of GW amplitude will ap- 
pear on the position at which the secondary BH moves to the po- 
sition of periastron. To some extent, we also find that the events 
that happened in maximum position are consistent with those of 
the optical outburst. For example, during 2005-2006 (Valtonen 
et al. 2008), the time of outburst from OJ287 is very close to 
a position of maximum GW radiation in our manifestation. In 
other words, if the OJ287 experiences an outburst in 2020, the 
very high possibility of GW detection can be investigated using 
pulsars around 2020. Thirdly, we find that the GW signal from 
this type of individual SMBH sources at the Earth is unlikely to 
be detected (the position of the symbol "today" is about 5 x 10~ 16 
in 2012.) using a current PTA detector, which agrees with esti- 
mates from previous studies (Hobbs et al. 2009, Sesana, Vecchio 
& Volonteri 2009; Yardley et al. 2010 ). Finally, these descrip- 
tions of the GW are the crucial criteria that we adopt throughout. 

2.2. Parameters in the Monte Carlo simulation 

To systematically investigate the GW radiation of galaxy pair 
systems, we performed a Monte Carlo simulation through which 
we followed the evolution of a sample of 1 million binaries. The 
description of our input physics parameters is given as follows. 
(i) The properties of the overall galaxy pair samples have been 
displayed in the released SDSS data (Ellison et al. 2010), includ- 
ing projected and angular separation, relative velocities (Av), 
stellar mass ratio, and redshift. We present polynomial curve fit- 
ting formulae, which we used to construct the distribution func- 
tion of galaxy pairs (see Figs. 1 and 2 of Ellison et al. 2010), 
which approximate the evolution of SMBHs for a wide range 
of mass M(10 8 - 1O 12 M ) and stellar environment log£(-1.4 - 
1.4Mpc~ 2 ). In our Fig. 3 (see panels (a)-(g)) we give the fitting 
curves of the physical properties for the galaxy pair sample us- 
ing these released SDSS data. 

(ii) The star formation rate(SFR) is taken to be constant in the 
simulation (Ellison et al. 2008). 

(iii) According to the derivation of the black hole mass function 
(BHMF) at low redshift (Tamura et al. 2006), a simple polyno- 
mial curve fitting approximation of BHMF is also used in panel 
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Fig. 2. Amplitude Aasa function of the observation time for the 
OJ287's model. The symbol "•" indicates the position at which 
the GW radiation from OJ287 is, ~ 5 x 10~ 16 in 2012. 



(h) of Fig. 3. 

(iv) We assumed that the merger rate of nearby (z < 0.1) galaxy 

pairs is 0.2 yr~' in the phase of coalescence with a chirp mass of 

10 10 M o . 
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Fig. 3. Fitting physical property curves of Monte Carlo parame- 
ters are derived from SDSS (Ellison et al. 2010) in panels (a)- 
(g), whose error bars represent the uncertainties of estimation 
due to the fitting errors in galaxy pairs of SDSS. The BH mass 
function curve is also display in panel (h), which is the same as 
that of Tamura et al. (2006). 



To construct the SMBH systems, other parameters used in 
the typical SMBH binary formation in this work are given as 
follows. 
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(i) The correlation between the mass Mbh of the BH and the this variation of the strain amplitude or luminosity with time 
velocity dispersion cr is given by (Tremaine et al. 2002) maybe provides an indirect evidence for a GW. Especially using 

20 years of observational time, the pulsar timing measurements 
log(MBH/M Q ) = a + /31og(cr/cr ), (14) will confirm a period variation in residual data. 

where cr = 200km/s, a = 8.13 and/3 = 4.02. 

(ii) For the BH mass and the luminosity relation, we adopted the 

B-band relation described by Bell et al. (2004) 



log(M BH ) = 1.19(logL- 10.0) + 8.18. 



(15) 



(iii) A gravitational potential in a galaxy includes the bulge, disk, 
and halo potentials. The disk and bulge potentials are derived 
from Miyamoto & Nagai (1975) 
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where R and z are the galaxy center coordinate parameters and 
2 The index i corresponds to bulge and disk. The 
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halo potential is given as (Paczynski 1990) 
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(17) 



A more detailed description of the parameters in Eqs. (16) and 
(17) can be found in Sect. 2 of Paczynski (1990). We obtained 
the eccentricity of the SMBH samples from the random num- 
ber (e = x, < x < 1). Throughout this paper, we adopted the 
cosmological model with Hq - 70kms~ 1 Mpc~ 1 , Qm = 0.3 and 
Q A = 0.7. 
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Fig. 4. Strain amplitude of GW radiation variation with time. 
The symbol "A" represents the average value position of the 
summed GW strain in each bin. The dotted line (log h - -15) 
marks the highest sensitivity of SKA. 



3. Results 

With a BPS simulation we obtained a series of SMBHB systems 
from a sample of 1 million initial SMBHBs. Based on the sensi- 
tivity curves of the PTA and SKA detectors (Yardley et al. 2010), 
we created the time series for the source signals and added the 
individual time series of the calculated sources to produce the 
total data stream. Similar calculations of the GW signal analysis 
combined with a GW detector can be found in Sesana, Vecchio 
& Colacino (2008), Liu (2009), and Ruiter et al. (2010). 

Figs. 4 and 5 show the variation of the strain amplitude and 
luminosity with time from the GW radiation of SMBHs dur- 
ing a long time observation (e.g. 20 years) in our Monte Carlo 
simulations. These tendency changes over time reveal that the 
GW signals obtained by current and planned detectors should 
be variable events. They will likely detect two higher GW ra- 
diation signals around 2020 and 2032, which is consistent with 
the predicted time of optical bursts of a merging galaxy (such as 
OJ287, Valtonen 2007). This is because we adopted a precession 
of orbital evolution of OJ287 as our orbital evolutionary model 
to investigate the orbital track changes caused by GW radiation 
in the simulations. In other words, according to our descriptions 
in Sect. 2, the maximum outbursts of GW radiation and the opti- 
cal outburst are happening at the exact time when the secondary 
BH crashes the periastron of the primary BH. 

Fig. 4 shows that the strain amplitude of GW radiation in 
the range of 4 x 10~ 17 to 7 x 10~ 14 will be detected by the SKA 
detector (Yardley et al. 2010). From Fig. 5 we can see that our 
calculated MBHB samples can radiate a high luminosity due to 
GWs of between ~ 1 X 10 48 and 10 51 erg-s~\ which possibly 
means that the influence of GW dominates the total energy loss 
of the galaxy pairs in the early phase of coalescence. Meanwhile, 
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Fig. 5. Similar to Fig. 4, but for the luminosity variation of GW 
radiations with time. 



Using a 20 year-integration time, we plot in Fig. 6 the GW 
radiation distribution from the resolved MBHB systems in our 
simulation. This shows that the predicted GW radiation is be- 
low the PTA sensitivity curve. Therefore the individual SMBH 
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Fig. 6. Strain amplitude Aas a function of the frequency for 
the BPS simulations of individual resolvable MBHBs. The gray 
scale gives the density distribution of the resolvable systems. 
Other lines are cited by the sensitivity studies (Sesana, Vecchio 
& Colacino 2008; Yardley et al. 2010) using the PTA detector. 

Table 1. The percentage of detectable sources with different ec- 
centricities (e). Note that the total number of detectable SMBH 
binaries in our model is 100. 



eccentricities 


< e < 0.2 


0.2 < e < 0.5 


0.5 < e < 0.7 


0.7 < e < 1.0 


percentages 


7.0% 


18.0% 


33.0% 


42.0% 



sources cannot radiate sufficient GW signals for the current de- 
tector. However, these signals will be detected by SKA, the 
detectable number of SMBHs is about 100, which is above 
the SKA detection threshold. Our prediction agrees well with 
that of Sesana, Vecchio & Colacino (2008), Sesana, Vecchio & 
Volonteri (2009), and Yardley et al. (2010) (see the discussion). 
Iwasawae al. (2010) have systematically investigated the im- 
portance of eccentricity in SMBH binaries due to GW radiation. 
For the 100 detectable sources in Fig. 6, we also give the per- 
centages of these sources for different eccentricities using Table 
1 . From it we find that the low-eccentricity binaries cannot dom- 
inate the detectable sources. The table also implies that once the 
individual SMBH binaries become highly eccentric, the GW ra- 
diation should be observable by the SKA detector. Meanwhile, 
we find for these detectable sources that the average of chirp 
mass (/Vlchiip) = ((m 1 m 2 ) 3/5 (m 1 +m 2 ) _1/5 > is about 8.9xlO ll) M . 

4. Discussion 

We have shown that there are good reasons to expect GW radi- 
ation from SMBHs and that some amplitudes will be measured 
by SKA in the range of 10~ 14 < h < 10~ 15 . The results from 
a BPS approach using Monte Carlo simulations are innovative 
in GW radiation studies, which also regulate the calculation to- 
gether with the analytical fitting functions of observation results 
from SDSS. In these close galaxy samples, very close SMBH 
binaries may produce very strong GW signals in orbital motion 
processes, which may be reflected by the profiles of GW forms. 
Especially the SMMB binaries with short orbital periods (<10 



yr) may even be detectable by photometric and high-energy radi- 
ation analysis using ground-based or space telescopes. Note that 
the precession of orbital period caused by GW radiation from 
highly eccentric SMBH binaries can influence the evolution of 
galaxy pairs. 

To study the SMBH systems such as GW sources, we needed 
some basic parameters for the SMBH systems (the primary BH 
mass ni2, the secondary BH mass mi, the orbital period P or b, and 
so on), which were derived from the distributions of SDSS data 
(Ellison et al. 2008, 2010) using the polynomial curve fitting ap- 
proximation. First, few studies focus on the GW radiation from 
SMBH systems together with the latest findings of observations, 
although our estimation methods are straightforward. Second, 
indeed, enough close galaxy pairs have been released from the 
SDSS data, but the error (or gray) bars (see Fig. 1 of Tamura et 
al. 2006 and Fig. 5 of Ellison et al. 2008) caused by the obser- 
vations still exist in our simulation. Therefore a maximum ratio 
with respect to observations is approximately 10, which is de- 
rived from the maximum difference obtained from error bars. 
Third, for the redshift evolution of SMBH systems, we inves- 
tigated the redshift distribution up to z ~ 0.1, which compares 
very well with a sky-averaged constraint on the merger rate of 
nearby SMBH sytems using a pulsar timing array measurement, 
implied by the low-redshift distance approximation r - cz/Hq. 

The main purpose of this study is to find out how many 
resolved SMBH sources we can explore with a pulsar timing 
detector from the GW radiation. Given the shape of the sky- 
and polarization-averaged sensitivity curve of the PTA (Yardley 
et al. 2010), the sensitivity curve is useful for comparing the 
strain amplitude of the GWs from an SMBH system with the 
average noise to see whether the individual sources can be de- 
tected by the PTA or SKA detector. A resolved source with fre- 
quency / and strain amplitude h that is observed over a time 
r o b S (=20 years) will appear in the Fourier spectrum of the data 
as a single spectral line. Therefore the SMBH systems, which 
are the only ones in their corresponding periodogram ranges 
(A/ = l/r o b s = 1.6 x 10~ 9 ) and have higher strains than the sen- 
sitivity noise curve of the detector, are called resolved sources. 

As we displayed in Fig. 6, the significantly different results 
can be explained as follows. Firstly, there are individual resolved 
SMBH systems in our nearby galaxy pairs above the predicted 
stochastic GW background, which is likely to be detected us- 
ing the PTA in the near future (Sesana, Vecchio & Colacino 
al. 2008). Note that our estimation of GW radiation agrees with 
that of Sesana, Vecchio & Colacino (2008), although these au- 
thors used different methods for the underlying SMBH popula- 
tion, for instance, we used a precession of orbital period evolu- 
tion of OJ287 and a series of derived observation results from 
SDSS in our simulation. Secondly, the distance was obtained us- 
ing D - cz/Hq, which is an acceptable approximation for rela- 
tively low redshift of galaxy pairs (Davis & Lineweaver 2004). 
All resolved sources in Fig. 6 are located at fairly low redshift 
(z < 0.1), and not at 0.2 < z < 1.5, which is one condition 
of Sesana, Vecchio & Volonteri (2009), who predicted that their 
individual SMBH is below the current sensitivity. Accordingly, 
the major difference in the estimating the resolved sources be- 
tween our model and the model of Sesana, Vecchio & Volonteri 
(2009) arises from the distance D. Thirdly, Iwasawa et al. (201 1) 
has pointed out that the eccentricity of most of SMBH bina- 
ries can result in a rapid merger through GW radiation. The 
75% detectable sources in Fig. 6 show high eccentricities, but 
Yarledy et al. (2010) did not consider the eccentric SMBH bina- 
ries. Consequently, the influence of the distance and eccentricity 
are the main difference from other papers. 
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If the orbital motion of the SMBH system is appropriately re- 
produced using PN2 methods, we expect that the maximum GW 
radiation should be associated with the precise evolution trace of 
each BH component movement in the full relativity field. This is 
because the secondary BH crashes the periastron, which leads to 
the outburst GW signal. Furthermore, based on our estimate of 
the merger rate of close galaxies and making several evolution- 
ary hypotheses, we can compute the number of SMBH sources 
that can be detected by the PTA (or SKA) detector. Additionally, 
the strength of GW signal from the merger of SMBH binaries 
will be blurred by other fierce energy release processes (such 
as supernova explosion, Mazzali et al. 2008). Finally, the evo- 
lutionary properties of SMBHs could indeed be more complex 
than we have described. For example, the mass transfer process 
between the components could be important and may even ac- 
count for quasars. Meanwhile, the merging process may occur 
more rapidly than the evolution of SMBH, and even more BHs 
may reside in the core of galaxy. 

To summarize, a numerical simulation of the GW radiation 
from the nearby galaxy pairs has been carried out using OJ287's 
modulation. We used the result of SDSS data and several notable 
conclusions from the literature to follow the SMBH evolution 
using a BPS approach. In addition to a GW detector (such as 
SKA), we used these methods to examine the GW variation with 
time and estimated the number of detectable SMBH sources us- 
ing observations made over 20 years. This study demonstrates 
the significant contribution to the GW radiation from the indi- 
vidual SMBH sources. 
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